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INTRODUCTION 

The deviation of a reflector antenna surface from a perfect 
parabolic shape causes degradation of the performance of the antenna. The 
shape of the antenna is therefore desired for several applications. If the 
shape of the antenna can be determined quickly during its manufacture, 
localized deviations from a perfect surface might be eliminated. If an 
antenna should become damaged, the location of the damage may allow 
easier repair. This is particularly important since the damage is not easy 
to see and is difficult to measure directly. 

The problem of determining the shape of the reflector surface in a 
reflector antenna using near field phase measurements in not a new one. A 
recent issue of the IEEE transactions on Antennas and Propagation (June, 
1988 ) contained numerous descriptions of the use of these measurements, 
including works by Y Rhamat-Samii, et al, W. Chujo, et al., and J. J. Lee, et 
al. These accounts use one of two methods: holographic reconstruction or 
inverse Fourier transform. 

Holographic reconstruction, used by Rahmat-Samii, makes use of 
measurement of the far field (amplitude and phase) of the reflector and 
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then applies the Fourier transform relationship between the far field and 
the current distribution on the reflector surface. 

Inverse Fourier transformation uses the phase measurements to 
determine the far field pattern using the method of Kerns. After the far 
field pattern is established, an inverse Fourier transform is used to 
determine the phases in a plane between the reflector surface and the 
plane in which the near field measurements were taken. 

These calculations are time consuming since they involve a 
relatively large number of operations. For the holographic reconstruction 
technique, the calculations are of the order of n 2 log(2n) floating point 
operations per phase measurement. The inverse Fourier transform method 
requires n 2 log(2n) calculations to obtain the far field pattern, followed by 
n 2 log(2n) operations to obtain the near filed phases again. 

A much faster method can be used to determine the position of the 
reflector. This method makes use of simple geometric optics to determine 
the path length of the ray from the feed to the reflector and from the 
reflector to the measurement point. This method takes only 57 floating 
point operations per phase measurement and gives the specular reflection 
point directly, rather than the phase at a plane near the reflector, as the 
inverse Fourier transform method does. 

For small physical objects and low frequencies, diffraction effects 
have a major effect on the error, and the algorithm provides incorrect 
results. It is believed (but not proven) that the effect is less noticeable 
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for large distortions such as antenna warping, and more noticeable for 
small, localized distortions such as bumps and depressions such as might 
be caused by impact damage. 

Determination of the applicable distortion feature sizes is outside 
the scope of this work. 


THE REFLECTOR SURFACE ESTIMATION ALGORITHM 
Necessary assumptions. 

The Reflector Surface Estimation (RSE) algorithm developed here, 
requires that there be no caustic points between the reflector surface and 
the measurement plane. If this assumption is met, each point on the phase 
measurement plane corresponds to either zero or one specular reflection 
point on the antenna surface. 

Geometry of the problem. 

The geometry used in the discussions throughout this document are 
shown in Figure 1. In accordance with normal conventions, the antenna 

radiates in the z-direction. A feed horn is located at the apparent focus 
(Xf,yf,Zf). A ray emitted from the feed intersects the reflector surface at 

the point (x,y,z). The ray is reflected from (x,y,z) and intersects the near 
field plane at a point (x a ,y a ,z z ). 
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Required data. 

The RSE algorithm requires transform phase measurements in the 
near field of a reflector antenna to the point on the antenna which caused 
the specular reflection. Required inputs to the basic algorithm are: 


Phase measurements 


Frequency 

Antenna feed location 


Reference length 
Phase measurement 


Absolute phase measurements or relative 
phase measurements which can be 
converted to absolute. 

The frequency at which the phase 
measurements were taken. 

The location of the antenna feed in the 
coordinate system in which the results are 
desired. 

One physical measurement which must be 
made to provide relative phase length. 

The distance from the origin of the 
coordinate system to the plane in which the 
phase measurements are made. 


Theory. 

The electrical distance from the feed of the antenna to the phase 
measurement plane can be found from: 


d — d re f - 


<kef - <t>ij 

k 


(1) 
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The distance consists of two components: the distance from the feed 
to the reflector and the distance from the reflector to the measurement 
point. 

dij = V (x-xf) 2 + (y-yf ) 2 + (z-z f ) 2 + V (x-Xa ) 2 + (y-y a ) 2 + (z-Za) 2 (2) 


It is desired to know the location of that reflector point. Since we 
know the phase at many points in the near field, we can calculate the 
angle of arrival of the ray. The partial derivatives are: 


_ foi+l,j ~ foi-l,j 
5x x i+l,j " x i-l,j 


dfoij _ foi.j+l ~ ^i.j-1 
9y yi.j+i - yi.j-i 


(3) 


m 


_ 1 d<t>ij 

x k 9x 

i 


m 


1 3<t>jj 
y k By t 


m, = - 


V 1 - m2 - 


(4) 


From these derivatives, we can define the path of the incoming ray: 


X = — ^-(z-Za) + X a 
m z ^ a 



(5) 


We define constants representing the slopes of the lines: 


Ci 




(6) 


Substituting (5) and (6) into (2) 

dij = V (Ci(z-Za) +X a -Xf ) 2 + (C 2 (z-Za) +y a -yf ) 2 + (z-Zf) 2 + V (x-Xa ) 2 + (y-ya ) 2 + (z-Za) 2 (7) 
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or: 


dij = V(-Ciz a +x a -x f + Ciz) 2 + (-C 2 z a +y a -yf + C 2 z) 2 + (z-z f ) 2 + V (x-xj 2 + (y-yj 2 + (z-zj 2 ^) 
Two new constants are now defined, 

di=-Ciz a + x a -Xf d 2 =-C 2 z a + y a -yf ( 9 ) 

Substituting into equation (8) 

dij = V(di + Ciz) 2 + (d 2 + C 2 z) 2 + (z-Zf) 2 + V (x-x a ) 2 + (y-yj 2 + (z-z a ) 2 (10) 

Expanding the first term and segregating powers of z, 

dij = V(d^ + 2diCiz + C?z2) + (d| + 2d 2 C 2 z + c£z 2 ) + ( z 2-2zz f + z?) + V(x-x a ) 2 + (y-y a ) 2 + (z-z a ) 2 

( 11 ) 


or: 


djj = V(d? + 4 + 4) + (2diCi + 2d 2 C 2 - 2z f )z + ( C? + C? + l)z 2 + V(x- Xa ) 2 + (y-yj 2 + (z-zj 2 

( 12 ) 


Three additional constants are defined: 


fi = d] + dj + zf 


( 13 ) 
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f2 = 2diCi + 2d2C2 - 2zf 


(14) 


f3 = C? + cl+l (15) 

Substituting (13) through (15) into (8) 

dij = Vfi + f 2 Z + f 3 Z 2 + V (x-Xa ) 2 + (y-ya) 2 + (z-Za) 2 (16) 

Substituting (6) into equation (16) yields 

dij = Vfi + f2Z + f 3 Z 2 + V Ci(z-Za ) 2 + C2(z-Za ) 2 + (z-Z a ) 2 (17) 

dij = Vfl + f 2 Z + f 3 Z 2 + V f 3 (z-Za ) 2 (18) 

Vfi + f2Z + f3Z 2 = dij - V f3(z-Za ) 2 ( 19 ) 

The second term on the right hand side of equation (19) can be either 
z-z a or z a -z. One root represents the desired solution and the other root 

represents a point along the ray but in the positive z direction from the 
near field plane. 

Squaring both sides and selecting the proper root, 
fi + f 2 z + f 3 z 2 = dfj - 2ff}( z-zjdjj + f 3 (z-z a ) 2 (20) 

Separating the powers of z, 

fi + f 2 Z + f 3 z 2 = dij - 2/TTzadij + f3Za + (2/FTdij - 2f 3 Za)z + f 3 Z 2 (21) 
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The 2 terms cancel, so 


fi + f 2 z = dfj - 2VT7z a dij + f 3 z} + (2VTJdij - 2f 3 z a )z 


( 22 ) 


Solving for z, 


fi - d-j + 2V7^z a djj - f 3 z a 

2/^dij - 2f 3 z a - f 2 (23) 


The x and y points may be found from equation (5). 


To obtain these results, the floating point operations in table 1 must 
be performed. 
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EFFECT OF NEAR FIELD GRID SIZE ON ACCURACY 


An attempt was made to determine the effect of the near field grid 
size on the accuracy. The accuracy should worsen with larger grid sizes 
because the partial derivatives are determined from the phases of the 
nearest neighbors, and, in the presence of distortion, the calculated 
partial derivatives differs from the true local partial derivative for large 
grid sizes. 


The analysis uses as a reflector model a parabola with cosine 
distortion in one of the axes. The surface is described by the equation 


z = 


x 2 + y 2 
2f 


+ 



Ymax “ y 
Ymax " Ymin 


with 8 ranging from 0 to 0.007 meters. 

The average error as a function of number of elements in the model 
antenna is plotted in figure 2 for several levels of distortion. As can be 
seen from the figure, the algorithm error varies nearly linearly with input 
distortion and is not greatly influenced by the element size. The 
invocations of RSE, and that the RSE algorithm failed for the large values 
of distortion for large grid sizes (evidenced by the curves which 
terminate early on the left of the plot). The RSE algorithm determined that 
some pairs of input phases were increasing or decreasing, wrote a 
message to the screen, and terminated the calculations for these cases. 
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For successful invocation of RSE, however, the accuracy of the result is 
not heavily influenced by grid size. 

This is not to say that the number of grid points is not an important 
parameter. If the number of grid points is small, the position of the 
reflector surface will be known at only a few points. 


EFFECT OF PHASE MEASUREMENT ERROR ON ACCURACY 

An important performance measurement for any algorithm that uses 
real measurements is the effect of errors in the measurements on the 
accuracy of the results. In order to determine the output of the program to 
input noise, Gaussian noise of various amplitudes was added to the input 
phase measurements. The results are shown in figure 3. 

The nonlinearity in the average output error as a function of input 
noise is apparently because the major effect on the error at low input 
noise levels is due to truncation errors in the algorithm. At higher levels, 
the error due to noise is the dominant part. 
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APPENDIX I 
AUXILIARY PROGRAMS 


This appendix contains description and listings of auxiliary 

programs used in the analysis. These programs include: 

vary: A program which uses all of the subroutines and functions 

below to produce error data based on variations in grid 
spacing, phase accuracy, and frequency. 

rnfgp: A subroutine which generates near field phase data on regular 

grid points based on user-supplied reflector distortion. 

fixphi: A subroutine which accepts the near field phase data supplied 

from rnfpg or from actual phase measurements and eliminates 
discontinuities which normally occur either at n and at -n or 
at 0 and 2ji. The input range is either (-n,n) or (0 .,tc) and the 
output range is unlimited. 

rseerr: A subroutine which includes the RSE algorithm and uses a 

user-supplied reflector distortion function (also supplied in 
rnfpg) to determine the error of the RSE algorithm. 

reffun: A function which is supplied to rnfpg and rseerr which returns 

the z position of a simulated reflector surface given the x and 


y coordinates. The partial derivatives with respect to x, y, and 
z are also given. 


Subroutine RNFPG. 

RNFGP theory. 

Program mfpg is an iterative procedure used to determine a point on 
a reflector, x r , y r> z r , which will reflect incident rays from a known feed 

point to a known point in a near field plane. Only two of the variables are 
required; the third can be determined because it is known that the point 
lies on the reflector surface. The projection of the geometry in the y=0 
plane is shown in figure 4. 

The procedure begins with the selection of a starting value for the 
solution. The assumption is made that the x and y coordinates on the 
reflector are close to the x and y coordinates in the near field plane. About 
this point, four rays are used to probe the location of the exact solution. 
These rays originate form the feed location, intersect the reflector at four 
points arranged about the assumed solution (see figure 5). 

Each of the rays is bounced off the reflector, following the laws of 
geometric optics. The intersection of the resulting ray and the near field 
plane is then calculated. These projections and the target grid point are 
shown in figure 5. 



From these points, a new value of x r , y r is selected by linear 
interpolation: 

x r ' - x r - 5 + 2 * 5 * (x a -x a3 ) / (x a1 - x a3 ) 
y r ' = y r - 8 + 2 * 8 * (y a -y a 4 .) / (y a 2 ‘ ya4) 


Up to this point, we have not discussed the selection of 5. Obviously, 
to converge to a solution, 8 must decrease with each iteration. The speed 
of convergence to a solution is directly related to the rate at which 5 
decreases. If 8 is decreased too quickly, however, x r ,y r may fall outside of 

the bundle of rays. This is usually not fatal, if the function is well 
behaved, but if it happens too often, divergence may occur. 

The parameter which will determine how quickly 8 can be reduced is 
related to the linearity of the mapping from the reflector position to the 
near field position. If the mapping is totally linear (e.g., no distortion), 
only one iteration is necessary. The more non-linearity, the more 
iterations will be needed. 

One rough indication of linearity can be obtained from the points 

already calculated. If the transformation were totally linear, the distance 
in the x axis from x a i to x a 2 would be the same as the distance in the x 
direction x a 4 to x a3 . Using the difference between the two distances 
divided by the total distance from x a i to x a3 as the measure of non- 
linearity, we have: 


skew x = |(x a1 + x a3 - x a2 - x a4 ) / (x a1 - x a3 )| 



A similar measure can be made in the y direction. 


skew y = |(y a1 + y a3 - y a2 - y a 4> 7 (Yal ' Ya3)l 

The program uses the non-linearity as the basis for the decrease in 
the spread of the packet of rays. The program starts with a 5 of 10% of 
the largest dimension of the antenna. After the first iteration, 5 is 
calculated from 


5 n+1 =5 n *(skew') 


Where 

skew' = max(0.1, min(0.9,max(skew x ,skewy)) 

The maximum value of skew x or skewy is used, so long as that value is 
greater than 0.1 and less than 0.9. 

After the new value for 5 is determined, a new bundle of rays is 
launched. A test is made to determine if the bundle of rays does enclose 

the solution. This can be determined by examining the intersection of the 
rays with the near field plane. x a3 should be less than x a and x a -| should 

be greater than x a , with similar requirements in the y axis. If one of these 

conditions is not met, an informative message is sent to the console and 
the value of 5 is automatically multiplied by 2. The iteration then 

continues. 
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After each iteration, a test is made to determine if the error in the 
near field plane has converged to within the maximum error used in the 
program's calling argument. If it has, the value is printed out to a file and 
the program continues with the next point in the near field plane. 


Subroutine fixphi. 

Subroutine fixphi takes as its input the results of rnfpg of near field 
phase measurements from an antenna facility and transforms the relative 
phase measurements ( -n < 0 < tt ) or ( 0 < 0 < 27 c ) to measurements which 
can be used to determine phase length. For example, if the following line 
were input into the program: 

1.0 1.5 2.5 0.5 1.5 

the program would convert the line to: 

1.0 1.5 2.5 3.64159 4.64159 

The program works by first examining the data to see if it meets one 
of the conditions: ( -k < 0 < % ) or ( 0 < 0 < 27 c ). If it meets neither 

condition, an error message is displayed on the console and the program 
terminates. If either condition is met, the program continues. 

The program continues by rewriting the input file and reading input 
while processing and printing the output. The first input value is special 
in that its value is always preserved. After the first value, each 
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measurement is examined to determine if it appears that the data has 
gone through a transition from -k ton or from 2k to 0. 

RNFPG performance. 

There were two figures of merit of the routine which were traded 
against each other to obtain maximum performance: computational speed 
and accuracy. Because of the iterative nature of the algorithm, additional 
accuracy can always be obtained by allowing more time for the 
computations, up to the precision limits of the machine. Double precision 
numbers were used as the default for the algorithm to limit the effect of 
machine precision on the output. 

Required accuracy is an argument in the invocation of RNFPG, and the 
algorithm will execute until that accuracy is obtained. 

In order to determine the effect of accuracy on expected execution 
time, number of iterations was plotted as a function of required accuracy. 
The results are shown in figure 6. 
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X 1 



Figure 1. Geometry used to develop the RSE algorithm 
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Tablel. Operations and timing. Times shown are for a Motorola 68881 co- 
processor operating at 40 mHz. Time is given in microseconds. 
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Figure 4. Problem geometry in the plane y=0 



Figure 5.rnfpg ray tracing 
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COMPUTER PROGRAM LISTINGS 


The programs used with this algorithm are listed on the following pages. 
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vary.f 


c vary.f compiles into a variety of programs depending on 

c the mode of compilation 

c Compilation must be done with the c preprocessor cpp. 

c One of the following may be defined 

c 

c (none) defaults to vary number of cells 

c ERROR varies the allowable error of rnfpg 

c 


program vary 

implicit double precision (a-h) 
implicit double precision (o-z) 
double precision xamp(8),yamp(8) 
common /partial/ pdfdx,pdfdy,pdfdz 
common /distort/ del , omega, xampl ,yampl 
integer type 

common /phys/ f ,xf ,yf , zf , zp,xmin,xmax,ymin,ymax, freq, type 
real maxerr(20, 20) , avgerr(20, 20) , rmserr(20,20) 
integer error(20,20) 
integer npts(8) ,nerrors(8) 

#def ine MAX_CURVE$ 8 
#def i ne HAX_POINTS 50 
#ifdef ERROR 

character*80 f i lename, plttt l ,xttl,yttl, zttl 
real xdata(50,8) 
real ydata(50,8) 
real avgitr 
common /perf/ avgitr 
#else ERROR 

real xdata(50,8) 
real adata(50,8) 
real rdata(50,8) 
real mdata(50,8) 
real sigma(50,8) 
real inacc(50,8) 

#endi f ERROR 

common /plot/ idist,igrid, 

1 maxerr(20,20) , avgerr<20, 20), rmserr(20,20),error(20,20) 
data xamp /0.00, 0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07/ 
data yamp /0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00/ 


fmaxerr=. 000001 

ymax=1 . 1 
ymin=0. 1 
xmax=0.5 
xmin=-0.5 
type = 1 

#ifdef ERROR 

inf i l el =5 1 3 

#else 

inf i l el =5 12 

inf i le2=( inf i le1*3)/4 

idiv = if ix( log(f loat( inf i lei ) )/log(2.0)-2.5) 
ngrids =2* if ix( log(f loat( inf i lei ) )/log(2.0)-2.5) 
inf i lei = inf i lel+1 
inf i le2=inf i le2+1 
#endif ERROR 


f = 1.0 


xf = 0.0 
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vary.f 


yf = 0.0 
zf = 1.0 

zp = 1.0 

periods =1.0 

omega = 2.0 * 3.14159265 * periods / (ymax*ymin) 

freq=305 00000000. 
c Wave number 

k = 2.0 * 3.141592650 / (300000000/freq) 


open(7,f i le="resul ts.dat") 

open(99,f i le= H contour.dat M ) 

minjdist=1 

max_dist=7 

do 10 i=min_dist,max_dist 
idist=i 
xampl=xamp( i ) 
yampl=yamp( i ) 
write(7, 105) 

105 formate 'avgerr.wpg' ) 
write(7,106) xampl,yampl 

106 formate 1 rnfpg error analysis 1 / 

1 1 amplitude of (x,y) distortion = e 1 

2 ,d18.10, 1 , 1 ,d18. 10, 1 > 1 ) 

x l ambda =0.1 

t = xampl 

del = t * xlambda 

#i fdef ERROR 

do 20 ierror=1 ,5 

ferrmax=.000001*ei0**e i error* 1 )) 

wri tee6,*)"max err = ",ferrmax 

wri tee7,*)"max err = ",ferrmax 

opene4, f i le="phl .dat", form=' UN FORMATTED 1 ) 

call rnfpge inf i lei , ferrmax) 

closee4) 

xdatae i error, i ) = ferrmax 
ydatae i error, i ) = avgitr 
20 continue 

nptsei) = ierror*1 

10 continue 

filename = "i ter2.wpg" 

pltttl = "RNFPG ITERATIONS VS ACCURACY" 

Xttl = "ACCURACY" 

yttl = "ITERATIONS" 

Zttl = "DISTORTION" 

call plotwpgef i lename,plttt l, xttl, yttl, zttl, 

1 . 000001 , . 01 , . 001 , 0 . , 10 . , 2 . , 1 , 0 , 

2 max_dist-min dist+1, 

2 npt$,xdata,ycfata) 

#else 

opene4, f i le="ph1 .dat", F0RM=‘ unformat ted* ) 

call rnfpge inf i lei , fmaxerr) 

closee4) 

openf3, f i le="ph1 .dat" , FORM=* unformatted* ) 

openf4,f i le="ph1 f .dat", F0RM=‘ unformatted 1 ) 

call fixphieinf i lei ) 

closee4) 

closed) 
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vary.f 


open(4, f i le="ph2 .dat" , FORM=' unformatted* ) 

call rnfpgCinf i le2,fmaxerr) 

close(4) 

open(3, f i le="ph2.dat", FORM =l unformatted' ) 

open(4,f i le="ph2f .dat l, ,FORH= , unformatted l ) 

call f ixphiCinf i le2) 

close(4) 

close(3) 

igrid=0 

do 20 j=1,idiv 

igrid=igrid+1 

xdata( Igrid, i ) = inf i le1/(2**( j * 1 )) 

open(4, f i le= 1 phi f .dat 1 , FORM= ' unformatted' ) 

call rsegrid( inf i lei , 1 + ( inf i lei - 1 )/(2**< j - 1 ))) 

adata( igrid, i )=avgerr( igrid, i ) 

mdata ( igrid, i ) =maxer r ( igrid, i ) 

rdata( igrid, i )=rmserr( igrid, i ) 

if(error( igrid, i) .ne.0)go to 21 

close(4) 

igrid=igrid+1 

xdata( igrid, i )=inf i Ie2/(2**(J* 1 )) 

open(4, f i le= 1 ph 2 f .dat 1 , FORM= 1 unformatted' ) 

call rsegrid( inf i le2,1+(inf I le2- 1 )/(2**( j-1))) 

adata( igrid, i )=avgerr( igrid, i ) 

mdata( igrid, i )=maxerr( igrid, i ) 

rdata( igrid, i )=rmserr( igrid, i ) 

if(error(igrid,i).ne.O)go to 21 

close<4) 

20 continue 

21 npts<i) = j*1 

do 30 j=1 , 10 

nerrors(j) = 10 

sigmaC j, i )=( (300000000. /freq)/(360.*. 01*5. )*f loat(j)) 
open(4, f i le=' phi f .dat 1 , FORM=' unformat ted* ) 
rewind(4) 

call rseerr(inf i le1,sigma(j', i), inacc( j, i ), ierrf Ig) 
close(4) 

30 continue 

10 continue 


call plotwpg("maxerr. wpg", 

1 "MAXIMUM ERROR AS A FUNCTION OF DISCRETIZATION", 

2 "NUMBER OF ELEMENTS ON A S IDE", "MAXIMUM ERROR", "DI STORTION" 

3 10., 1000., 10., 0.,. 00003, .00001,1,0, 

4 maxjdist-min dist+1, 

5 npts,xdata, mdata) 

call plotwpg("avgerr.wpg", 

1 "AVERAGE ERROR AS A FUNCTION OF DISCRETIZATION", 

2 "NUMBER OF ELEMENTS ON A SIDE", "AVERAGE ERROR", "DISTORTION" 

3 10. ,1000., 10., 0.,. 00003,. 00001, 1,0, 

4 max_dist-min dist+1, 

5 npts,xdata, adata) 

call plotwpgO'rmserr.wpg", 

1 "ROOT MEAN SQUARE ERROR AS A FUNCTION OF DISCRETIZATION", 

2 "NUMBER OF ELEMENTS ON A SIDE", "RMS ERROR", "DISTORTION", 

1 10., 1000., 10., 0.,. 00003,. 00001, 1,0, 

2 max_dist*min dist+1, 

3 npts,xdata,rdata) 

call plotwpg( n randerr.wpg n . 
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vary.f 


1 "RSE AVG ERROR IN THE PRESENSE OF NOISE" 

3 oTolKffi3*o®:^6^ ERR0R ' HETERS,, ' ml * 

4 max_dfst-mfn_dist+1 , 

5 nerrors, sigma, inacc) 

#endif ERROR 
stop 
end 
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reffun.f 


c Calculate the z coordinate of the reflector surface 

real*8 function reffun(x,y) 

implicit real*8 <a-h) 
implicit real*8 (o-z) 
real argx, sinfunx, cosfunx 
real argy, sinfuny, cosfuny 
integer type 

real*8 x,y, f, del , omega ,ymax 
real*8 temp 

common /partial/ pdfdx,pdfdy,pdfdz 
common /distort/ del, omega 

common /phys/ f ,xf ,yf # zf ,zp,xmin # xmax,ymin,ymax,freq, type 

c x and y are the x and y positions of the point 
c del is the distortion amplitude factor 
c omega is the distortion wave number 
c ymax is the maximum y value 


if (type . eq. 1) then 

c This type is for an antenna with sinusoidal distortion that varies 
c only with the y variable 

argy = omega * (ymax-y) 
sinfuny = sin(argy) 
cosfuny = cos(argy) 
temp = -0.5/f 
pdfdx “ x*temp 

pdfdy = y*temp - del * omega * sinfuny 
pdfdz a 1,0 

reffun = (x**2+y**2)/(4.0*f )+del*cosfuny 
return 

else if (type .eq. 2 ) then 

c This type is for an antenna with sinusoidal distortion that varies 
c only with the x variable 

argx = omega * (xmax*x) 
sinfunx = sin(argx) 
cosfunx = cos(argx) 
temp = -0.5/f 

pdfdx = x*temp - del * omega * sinfunx 
pdfdy = y*temp 
pdfdz =1.0 

reffun = (x**2+y**2)/(4.0*f )+del*cosfunx 
return 

else if (type .eq. 3 ) then 

This type is for an antenna with sinusoidal distortion that varies 
with both the x and the y variable 

argx = omega * (xmax-x) 
argy = omega * (ymax-y) 
sinfunx = sin(argx) 
cosfunx = cos(argx) 
temp = -0.5/f 

pdfdx = x*temp - del * omega * sinfunx 
pdfdy = y*temp - del * omega * sinfuny 
pdfdz =1.0 

reffun = (x**2+y**2)/(4 . 0*f )+del*cosfunx +del*cosfuny 
return 
end if 
end 
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subroutine rnfpg(nphasegp, ferrmax) 

c The purpose of this program is to detect an antenna 
c reflector surface from the near field phase distribut ion. 
c The near field phase distribution is defined on an nxn 
c rectangular grid system. 

c Feeder location(xf # yf ,zf ) , value of lambda, diameter of 
c the reflector aperture(d) 

c ******************* *************** ***************************** 

implicit real*8 (a-h) 
implicit real*8 (o-z) 
real*8 xa, ya 
real avgitr 

common /partial/ pdfdx,pdfdy,pdfdz 
integer type 

common /phys/ f,xf ,yf ,zf , zp,xmin, xmax,ymin,ymax, freq, type 
common /per f/ avgitr 
dimension phi (1026) 
open(11,f i le^'res* ) 

pi = 3.1415926 

c Diameter of the region of interest on the reflector 
d=1 . 

c Near field grid spacing 
delnf = d/(nphasegp- 1 ) 

k = 2*3.1415926 /(300000000. / freq) 

fctrmin = 0.1 
fctrmax - 0.9 

nf i rst=1 

wri te(6, 998) ferrmax 
998 formatO’max err =",g8.3) 

nloops=0 

c Scan through the x axis of the near field 
do 10 i=1,nphasegp 

c call tick(i) 

c Scan through the y axis of the near field 

do 20 j=1,nphasegp 

c Define the x,y,z coordinates in the near field 

xa = xmin + delnf * (i-1) 
ya = ymin + delnf * (j*1) 
za = zp 

c set up for search in reflector plane 

xr=xa 
yr=ya 

delta=dmax1(dabs(xmax-xmin),dabs(ymax*ymin))/10. 

c reflector plane iteration loop 

c 

601 continue 

nloops=nloops+1 

c 

c launch four rays 

call qray(xr+del ta,yr,xa1 ,ya1) 
call qray(xr,yr+delta,xa2,ya2) 
call qrayCxr-del ta,yr,xa3,ya3) 
call qray(xr,yr*delta,xa4,ya4) 
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649 
648 

650 

1 

1 

652 


if Cxa3.gt.xa)goto 649 
ifCxa.gt.xal )goto 649 
if Cya4.gt.ya)goto 649 
ifCya.gt.ya2)goto 649 
go to 650 

cont i nue 
wri te(6,648) 

formate got outside of bundle of rays’) 
del ta=delta*2 
goto 601 

continue 
if(xa1 ,ne.xa3) 

xr=2*del ta/(xa1 -xa3)*(xa*xa3) + xr * delta 
if Cya2.ne.ya4) 

yr=2*delta/Cya2-ya4)*Cya-ya4) + yr * delta 
conti nue 


skewx = 0 
skewy - 0 

ifCCxal • xa3) .eq. 0)goto 680 

skewx = dabsCC(xa1+xa3)*Cxa2+xa4))/Cxa1 -xa3)) 
i f C(ya2 * ya4 ) .eq. 0)goto 681 

skewy = dabsCC(ya2+ya4)*Cya1+ya3))/Cya2*ya4)) 
factor=2*dmax1 Cskewx, skewy) 
factor=dmax1 (fctrmin, factor) 
factor=dmin1 ( fctrmax, factor) 

del ta=del t a* factor 

if CCabsCCxal -xa3)*pdfdxf )+absCCya2-ya4)*pdfdy)) 

1 .gt.ferrmax) goto 601 

zr=ref funCxr,yr) 

di st=d$qrt C Cxf *xr)**2+Cyf -yr)**2+Czf - zr)**2)+ 

1 dsqrtC (xa-xr)**2+Cya-yr)**2+(za-zr)**2) 

phi ( j)=k*di st 
nphi=ph i C j )/(2*pi ) 
ph i C J ) =ph i C J* ) * nph i *2*pi 
i f (nf i rst .ne. 1 )goto 200 
write(4)dist,phi( j ) 

199 format(f14.8/f 14.8) 
nf i rst=0 

200 continue 

20 continue 

wri teC4)Cphi ( j ) , j=1 , nphasegp) 

555 formatCl026f 14.8) 

10 continue 

close(ll) 

avgitr = f loatCnloops)/f loat(nphasegp**2> 
write(7,997)avgitr 
wri te(6, 997)avgi tr 
997 formatC'avg iter = ”f8.3) 
return 
end 


680 

681 


subroutine qray(xr,yr,xa,ya) 
implicit real*8 Ca-h) 
implicit real*8 (o-z) 
common /partial/ pdfdx,pdfdy,pdfdz 
integer type 

common /phys/ f ,xf ,yf , zf , zp^xminjxmax^ymin^ymax, freq, type 
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real*8 Ux, I2x, Uy, I2y, 1 1 z , I2z 
zr = reffun(xr,yr) 

L lx = xf - xr 
Uy = yf - yr 
Uz = zf • zr 

absnrml = pdfdx**2 + pdfdy**2 + pdfdz**2 

r = 2.0*(pdfdx * llx + pdfdy * Uy + pdfdz * Uz)/ab$nrml 

I2x = Ux - r * pdfdx 

I2y = Uy • r * pdfdy 

I2z = Uz • r * pdfdz 

quick = (zp-zr)/l2z 

xa = quick * I2x + xr 

ya = quick * I2y + yr 

return 

end 
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subroutine f ixphi(ninf i le) 
implicit double precision (a*h) 
implicit double precision (o*z) 
double precision phi (1026) 

pi = 3.14159265 
twopi = 2*3.14159265 

read(3)dref ,phi ref 
wri te(4)dref ,phi ref 
199 format(f 14.8/f 14.8) 
i first = 1 
irow = 0 

do 10 i=1 ,ninf i le 

read(3)(phi (n),n=1 # ninf i le) 

if ( if i rst .eq. 1) phi row = phi(1) 

i first = 0 

phi last = phi row 

if((phi(l)*phirow) .It. (-1*pi)) then 
irow - irow + 1 
phi last=phi (1 ) 
phi row=phi ( 1 ) 

else if((phi(1)-phirow) .gt. pi) then 
irow » irow * 1 
phi last=phi(1) 
ph i row=phi (1 ) 

endi f 

icol = irow 
do 20 j=1 , ninf i le 

if ((phi( j) -phi last) .It. (*1*pi)) then 
icol = icol + 1 

else if((phi( j)-philast) .gt. pi) then 
icol = icol - 1 

endi f 

phi last=phi( j) 

phi(j) = phi(j) + icol*twopi 
20 continue 

wri te(4)(phi (n),n=1 ,ninf i le) 

10 continue 

7 format(1026f 14.0) 

8 format(1026f 14.8) 
return 

end 
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subroutine rsegrid(ninfi le, ntouse) 

c the purpose of this program is to detect an antenna reflector 

c surface from the near field phase distribution. The near field 

c phase distribution is defined on an nxn rectangular grid system 
c 

c surface detection 

implicit double precision (a-h) 
implicit double precision (o*z) 
double precision lambda,mx,my,mz,k 
double precision phi(3,1026) 
double precision dummy(1026) 
real*8 reffun 
integer contour( 1026) 
common /partial/ pdfdx, pdfdy,pdfdz 
common /distort/ del , omega, xampl ,yampl 
integer type 

common /phys/ f ,xf ,yf , zf ,zp,xmin,xmax,ymin,ymax, f req, type 
real maxerr(20,20) ,avgerr(20,20) ,rmserr(20,20) 
integer error(20,20) 
common /plot/ idist,igrid, 

1 maxerr(20,20),avgerr(20,20),rmserr(20,20) f error(20,20) 

character*4 comment 


error< igrid, idist)=0 

comment = 1 5 

errmax = 0. 

sumer r = 0. 

sumsq = 0. 

nsum = 0 


i f (ntouse. It. 100) write(99,801 )ntouse # idist 

801 formate e‘ , i4.4, i2.2 , ' .dat 1 ) 

if (ntouse, It . 100) wr i te(99,802)ntouse,del 

802 formate 1 ALGORITHM ERROR WITH ',14.4, 1 GRIDS AND ', 

1 F9.3, 1 MM MAX SIN DISTORTION') 

i f (ntouse. I t . 100) wri te( 99,803) 1 , 'err<. 00003* 
i f (ntouse. I t . 100) wri te( 99, 803)2, *err<. 00002' 
if (ntouse. It. 100) wri te(99, 803)3, ' err<. 00001 1 
if (ntouse. It .100) wri te( 99, 803)4, ' err<. 00000' 
if(ntouse.lt. 100) wri te(99, 803)5, 'err>. 00000' 
if(ntouse. It. 100) wri te(99, 803)6, 'err>. 00001 1 
if (ntouse. It. 100) wr i te(99, 803)7, 'err>. 00002' 
if (ntouse. It. 100) wri te( 99, 803)8, 'err>. 00003' 

803 formate \i3,' \a) 

if (ntouse. It. 100) wri te(99,804)ntouse-2 

804 formate 1 , i 3 ) 


za = zp 

pi = 3.14159265 

read(4)dref , phi ref 

lambda=300000000./f req 
k = (2*pi )/ lambda 
delx=(xmax-xmin)/(ntouse- 1 ) 
dely=(ymax*ymin)/(ntouse- 1 ) 


ntoskip=(ninf i l e/( ntouse* 1 ) )*1 
write(6,*)' ntoskip = ',ntoskip 


read in two lines of input to start the process 
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if(ntoskip.eq.O) then 
read(4)(phi(2,n),n=1 .ntouse) 
else 

read(4)phi(2, 1), 

1 ((dummy(iskip), iskip=1 , (ntoskip)),phi(2 # n),n=2 # ntouse) 
endi f 

if (ntoskip.ne.O) then 
do 880 iskip=1 . (ntoskip) 

880 read(4)dunmy(1 ) 
endi f 

i f (ntoskip. eq.O) then 
read(4) (phi(3,n),n=1 .ntouse) 
else 

read(4)phi(3,1), 

1 ((dummy(iskip), i ski p=1 , (ntoskip)), phi (3, n),n=2, ntouse) 
endi f 

i f (ntoskip.ne.O) then 
do 881 iskip=1 .ntoskip 
881 read(4)dummy( 1 ) 

endi f 

do 400 n=1, ntouse 
phi (2.n)=phi (2. n) -phi ref 
400 phi(3,n)=phi(3.n) -phi ref 
do 10 i=2,ntouse-1 
c wri te(7.*)sumerr 
xa=xmin+( i - 1 )*delx 

c prepare to read in a new row 

do 11 i 1=1. ntouse 

phi ( 1 , i 1 )=phf (2, i 1 ) 

11 phi (2. i 1 )=phi(3, il ) 


i f (ntoskip. eq.O) then 

read(4)(phi (3,n),n=1 .ntouse) 

else 

read(4)phi(3,1), 

1 ((dummy(iskip). i ski p=1, (ntoskip)), phi (3, n),n=2. ntouse) 
endif 


if (ntoskip. ne.O.and. i . ne. (ntouse- 1)) then 
do 882 iskip=1 .ntoskip 
882 read(4) dummy (1 ) 

endi f 

do 401 n=1 .ntouse 
401 phi(3,n)=phi(3,n)-phiref 

do 20 j=2, ntouse- 1 
ya=ymin+( J * 1 )*dely 
d=dref+1/k*(phi (2, j )) 

c dref and phi ref are measured quantities 

dph i dx=ph i ( 3 , j ) - ph i ( 1 . j ) 
i f ( dph i dx . g t . pi ) dph i dx=dph i dx - 2*pi 
if(dphidx. It . ( - 1*pi ))dphidx=dphidx+2*pi 
i f ((dphidx. It . - 1 ) .or. (dphidx.gt.1 ))then 
if (error ( igrid, idist) .eq.O) then 
error ( igrid, idist)=1 
comment=' ?? 1 
wri te(6.570) 
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endi f 

end if 

570 formate grid size too Large') 

dph i dx=dph i dx/ ( 2*de lx) 

dph i dy=ph i ( 2 , j* + 1 ) * ph i ( 2 f J - 1 ) 
i f ( dph i dy . g t . p i ) dph i dy=dph i dy • 2*pi 
i f (dph idy . 1 1 . ( - 1 *p i ) )dph i dy=dph i dy+2*pi 
i f ((dphidy. It . * 1 ) .or. (dphidy. gt. 1))then 
if(error( igrid, idist) .eq.0)then 
error ( l grid, idist )=1 
comment = ' ?? 1 
wri te(6,570) 
endi f 

end i f 

dph i dy=dph i dy/ ( 2*de L y ) 
mx = dphidx/k 

my = dphidy/k 

mz = dsqrt ( 1 . 0-mx**2-my**2) 
cl = mx/mz 

c2 = my/mz 

dl = -c1*za + xa - xf 

d2 = -c2*za + ya - yf 

fl = d1**2 + d2**2 + zf**2 
f2 = 2.0 * (dl*c1 + d2*c2 - zf) 
f3 = c1**2 + c2**2 + 1.0 

znu = -fl + (f3*za**2 - 2.0 * dsqrt(f3)*d*za 
1 + d**2) 

znd = < f2 + f3*2.0*za * 2.0 * dsqrt(f3)*d) 
z=znu/znd 
x=c1*(z-za) + xa 
y=c2*(z’za) + ya 
c wri te(6,991 )x,y,z 

991 formate x= ”, f 18. 10, ", y= ", f 18.10, ", z= ",f18.10) 

err = z - reffun(x,y) 
if (err.gt .0)then 

contourCJ )=5 

i f (err.gt . .0001 )contour( ] )=6 
if (err.gt. .00Q2)contour( j)=7 
l f (err.gt . .0003)contour( j )=8 

else 

contour( j)=4 

if (err. It. * .0001 )contour( j )=3 
if (err. It. ■ ,0002)contour( J )=2 
i f (err. It.* .0003)contour( j )=1 

endi f 

derror= dabs(err) 

c write(6,556)x,y,z,derror 

556 format(5f 14.8) 

if(derror .gt. errmax) errmax=derror 
sumerr = sumerr + derror 
sumsq = sumsq + derror**2 
nsifm = nsum + 1 
20 continue 

if (ntouse. I t . 100) wri te(99,819)(contour< j ) , J=2, ntouse* 1 ) 
819 formate 1026i 1 ); 

10 continue 

rms = dsqrt (sumsq/nsum) 
average = sumerr/nsum 
maxerr( igrid, idist)=errmax 
avgerrC igrid, idist)=average 
rmserr( igrid, idist)=rms 
wri te(7, 557) errmax, average, rms 

557 formate maxerr = 1 , d14.8, 

1 1 sumerr/n = d14.8, 

2 ' sumsq/n = d14.8) 

wri te(7,558)ntouse, average 

558 formate ',14,' *,f14.8) 

return 

end 
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subroutine rseerr(ninf i l e, sigma, inacc, error) 
c 
c 

c surface detection 

implicit double precision (a-h) 
implicit double precision (o-z) 
real sigma, inacc, gasdev 
integer error 

double precision lambda, mx, my, mz,k 
double precision phi (3, 1026) 
real*8 reffun 

common /partial/ pdfdx, pdfdy,pdfdz 
cormon /distort/ del , omega, xampl ,yampl 
integer type 

common /phys/ f ,xf,yf , zf ,zp,xmin,xmax,ymin,ymax,freq,type 
logical exist, opened 
integer ios,nr 

integer argl ,arg2, arg3,arg4,arg5 


1 continue 
error=0 
errmax = 0. 
sumerr = 0. 
sumsq = 0. 

nsum = 0 

za = zp 


pi = 3.14159265 

read(4,end=699,err=698)dref ,phi ref 

lambda=300000000./f req 
k = (2*pi )/lambda 
delx=(xmax-xmin)/(ninf i le- 1 ) 
dely=(ymax-ymin)/Cninf i le* 1) 
c 

c read in two lines of input to start the process 
c 

read(4,end=699,err=698)(phi (2,n) ,n=1 ,ninf i le) 
read(4,end=699,err=698)(phi (3,n) ,n=1 ,ninf i le) 

do 400 n=1 , ninf i le 

perturb = sigma * gasdevC) 
phi<2,n)=phi(2,n) -phi ref + perturb 
perturb = sigma * gasdevC) 
phi (3,n)=phi (3,n) -phi ref + perturb 
400 continue 

do 10 i =2, ninf i le-1 
xa=xmin+( i -1 )*delx 

c prepare to read in a new row 

do 11 i 1=1, ninf ile 

ph i ( 1 , i 1 ) =ph i ( 2 , i 1 ) 

11 phi (2, il )=phi(3,i1) 

read(4,end=699,err=698)(phi(3,n),n=1 ,ninf i le) 

do 401 n=1 ,ninf i le 
perturb = sigma * gasdevC) 
phi(3,n)=phi (3,n) -phi ref + perturb 
401 continue 

do 20 j=2,ninf i le- 1 
ya=ymin+( j-1)*dely 
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d=dref+1/k*(phi (2, j )) 

c dref and phi ref are measured quantities 

dph idx=ph i (3, j ) -ph i ( 1 , j ) 
if(dphidx.gt .pi )dphidx=dphidx-2*pi 
if(dphidx. It . ( * 1*pi ))dphidx=dphidx+2*pi 
if((dphidx. Lt. - 1 ).or.(dphidx.gt.1))then 
if (error. eq.O) then 
error=1 
write(6,510) 

510 formate 'grid size too large 1 ) 

endi f 

endi f 

dph i dx=dph i dx/ ( 2*de l x ) 

dphidy=phi(2, j+1)-phi(2, j*1 ) 
if(dphidy.gt.pi )dphidy=dphidy-2*pi 
i f (dphidy. I t . ( * 1*pi ))dphidy=dphidy+2*pi 
i f( (dphidy. It . *1 ) . or. (dphidy. gt. 1 ))then 
if (error. eq.0)then 
error=1 
wri te(6,510) 
endi f 

end i f 

dph i dy=dph i dy/ ( 2*de l y ) 
mx = dphidx/k 

my = dphidy/k 

mz = dsqrt ( 1 . 0-mx**2-my**2) 
cl = mx/mz 

c2 = my/mz 

dl = -c1*za + xa - xf 

d2 = -c2*za + ya - yf 

fl = d1**2 + d2**2 + zf**2 
f 2 = 2.0 * (d1*c1 + d2*c2 • zf) 
f3 = c1**2 + c2**2 + 1.0 

znu = -fl + (f3*za**2 - 2.0 * dsqrt(f3)*d*za 
1 + d**2) 

znd = (f2 + f3*2.0*za - 2.0 * d$qrt(f3)*d) 

z=znu/znd 

x=c1*(z*za) + xa 

y=c2*(z-za) + ya 

err = z - reffun(x,y) 
derror= dabs(err) 

if(derror .gt. errmax) errmax=derror 
sumerr = sumerr + derror 
sumsq = sumsq + derror**2 
ns urn = nsum + 1 
20 continue 

819 format(1026i 1 ); 

1 0 cont i nue 

rms = dsqrt(sumsq/nsum) 
average = sumerr/nsum 
inacc = average 

wri te(7, 555)sigma, inacc, errmax, average, rms 
555 formate , sigma= l ,g!4.8, 1 , inacc=' , g14.8, 1 ^rrmaxs 1 , 914.8, 

1 averages', d14.8,',rms= d14.8) 

write(6,557)sigma, inacc 

557 formate 'sigmas' , g14.8, 1 , inacc=‘ ,d14.8) 
write(7,558)ninf i l e, average 

558 formate *,i4,' 1 ,f 14.8) 

return 

698 write(6,*) ,, -**error---" 

return 

699 write(6,*)"- - -end- * 

call ERRSNS(arg1 , arg2,arg3,arg4,arg5) 
write(6,*)arg1 ,arg2,arg3,arg4,arg5 

inqui re(4,EXI$T=exist,NEXTREC=nr, IOSTAT=ios,OPENED=opened) 

write(6,*)"exist", exist ,"ios tat", ios^opened' 1 , opened 

wri te(6,*)"next record = ",nr 

wri te(6,*)"n=",n, ,, ,nsum= ,, ,nsum, ,l f i=" i 

call flush(6) 

rewind (4) 

goto 1 
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program vet r 

c this program reads the coordinates of any point on the reflector 

c and determines the cooresponding coordinates on the near field 

c plane 

implicit real*8 (a*h) 
implicit real*8 (o-z) 

dimension x(9,9) ,y(9,9), 2(9,9) ,k(81 >, 1(81 ),xa(9,9),ya(9,9) 

real*8 l lx, lly, llz, l 1 ,ul1x,ul1y,ul 1 z 

real*8 r,ul2x,ul2y,ul2z 

real*8 pdfdx,pdfdy,pdfdz, absnrml 

real*8 nrmlx,nrmly, nrmlz 

open(12,f i l e^ 1 data 1 > 

open(13, f i l e= 1 res 1 , status= , old‘ ) 

write(6,*) ‘enter distortion factor t 1 

read(5,*) t 

x lambda = 0.1 

nperiod = 1 

f = 1.0 

ymax = 1.1 

ymin =0.1 

d = 1. 

delref = 0.125 
zp = 1. 

del = t * x lambda 

xf = 0.0 

yf = 0.0 

zf = 1.0 

pi = 3.14159265 

omega = 2.0 * pi * nperiod/ (ymax* ymin) 
c= 2 * pi /x lambda 
zaa = zp 

nrefgp = d/del ref + 1 
do 10 i=1,nrefgp 

do 20 j=1, nrefgp 

read(13,3)k(i), l(i),x(i, j),xa(i , j),y(i, j),ya(i, j),z(i, j) 
3 format(2i3,5f 11 .7) 

l lx = xf * x(i , J) 
lly = yf * y(i,i) 
llz = zf * z(i,J) 

11 = sqrt ( l 1x**2 + I1y**2 + I1z**2) 
pdfdx = -x(i, j)/(2.0*f) 

pdfdy = -y(i, j)/<2.0*f ) * de l *omega*sin( omega* (ymax 

1 -V(iJ))) 

pdfdz =1.0 

absnrml = sqrt(pdfdx**2 + pdfdy**2 + pdfdz**2) 

nrmlx = pdfdx/absnrml 

nrmly = pdfdy/absnrml 

nrmlz = pdfdz/absnrml 

ul lx = Ux / 11 

ully = lly / 11 

ullz = llz / 11 

r = 2.0*(nrmlx*ul1x+nrmly*ul1y+nrmlz*ul1z) 
ul2x = ullx - r*nrmlx 
ul2y = ully - r*nrmly 
ul2z = ullz - r*nrmlz 
c in original z(i,j) was z 

xaa = (zaa-z( i , j))*(ul2x/ul2z) + x C i , j > 
yaa = (zaa- z( i , j) )*(ul2y/ul2z) + y(i,j) 
wr i te( 1 2 , 5 )k< 1 ) , l ( j ) , x( i , j ) , xaa ,y( i , j ) , yaa , z( i , j ) , zaa 
5 format(2i3, 6f 10.5) 

20 continue 

10 continue 

close(13) 
close(12) 
end 
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function gasdevO 

c returns a normally distributed deviate with zero mean 

c and unit variance 

double precision rand 
data i set/0/ 
data gset/0.0/ 
if( iset .eq.O) then 
1 vl = 2. * rand() - 1. 
v2 = 2. * rand() - 1. 
r = v1**2 + v2**2 
if(r.ge.1)go to 1 
fac = sqrt(-2.*log(r)/r) 
gset = vl * fac 
gasdev = v2 * fac 
iset = 1 
else 

gasdev = gset 
iset = 0 
endi f 
return 
end 
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program vary 

implicit double precision (a-h) 
implicit double precision (o-z) 
double precision xamp<8) / yamp(8) 
common /partial/ pdfdx,pdfdy,pdfdz 
common /distort/ del , omega ,xampl f yampl 
integer type 

common /phys/ f ,xf ,yf # zf , zp,xmin,xmax,ymin,ymax, freq,type 
real maxerr(20,20) ,avgerr(20,20) , rmserr(20,20) 
integer error(20,20) 

integer npts(8) ,nerrors(8) 

#def i ne MAX CURVES 8 
^define MAX_POINTS 50 

c character*80 f i lename,pltttl,xttl,yttl,zttl 

tfifdef ERROR 

real xdata(50,8) 
real ydata(50,8) 
real avgitr 
common /perf/ avgitr 
#else ERROR 

real xdata(50,8) 
real adata(50,8) 
real rdata(50,8) 
real mdata(50,8) 
real sigma(50,8) 
real inacc(50,8) 

#endif ERROR 

common /plot/ idist,igrid, 

1 maxerr(20, 20) ,avgerr(20,20), rmserr(20, 20),error(20,20) 

data xamp /0.00, 0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07/ 
B 

data yamp /0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00/ 

integer at ime(3) , iday, i month, iyear 
500 formate i 2. 2, V 1 , i2.2, ‘/' , i2.2, ' ', i2.2, ' : « , 1*2.2* . ' i2.2) 

call itime(atime) 
call idate( iday, i month, iyear) 
wri te(6,500)imonth, iday, iyear, a time 
B 

fmaxerr=. 000001 

ymax=1 . 1 
ymin=0 . 1 
xmax=0.5 
xmin=*0.5 
type - 1 

inf i l el =1024 

idiv = if ix(log(f loatCfnf f le1))/log(2.0)-2.5) 
ngrids =2* if ix(log(f loat(inf i lei ))/log(2.0)*2.5) 
inf i le1=inf i lei +1 


f = 1.0 


xf = 0.0 
yf = 0.0 
zf = 1.0 

zp = 1.0 

periods = 1 .0 

omega = 2.0 * 3.14159265 * periods / (ymax-ymin) 

freq=30500000000. 
c Wave number 
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vary 1 k.f 


k = 2.0 * 3.141592650 / (300000000/f req) 


open(7, f i le=" results.dat") 

open(99,f ile= n contour .dat") 

min_dist=1 
max dist=7 


read(5,997) i 

997 formate i3) 

if(i . lt.min_dist.or. i .gt.max_dist) then 
write(6,999) i 

999 formate* bad distortion selector: 1 *, i4); 

stop 

else 

wri te(6,998) i , xamp(i), yamp(i) 

998 format ("doing i teration #'*, i2,"; ",g10.4,»,»,g10.4) 
endif 


c do 10 i =min_dist,max_dist 

idist=i 
xampl=xamp( i ) 
yampl=yamp( i ) 
write(7, 105) 

105 formate 'avgerr.wpg* ) 
write(7,106) xampl,yampl 

106 formate rnfpg error analysis 1 / 

1 1 amplitude of (x,y) distortion = (' 

2 ,d18. 10, « , 1 ,d18. 10 # 1 ) « ) 

x lambda = 0.1 

t = xampl 

del * t * x lambda 

open(4, f i le="ph1 .dat" , FORM= 'unformatted* ) 

call rnfpgC inf i lei , fmaxerr) 

close(4) 

open(3, f i l e="ph 1 .dat", FORM=’ unformat ted' ) 

open(4, f i l e="ph1f .dat", F0RM=* unformat ted* ) 

call f ixphi(inf i le2) 

close(4) 

close(3) 

21 npts(i) = j*1 

10 continue 

call itime(atime) 
cal l idate( iday, imonth, iyear) 
write(6,500)imonth, iday, iyear, atime 
stop 
end 
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